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Abstract 

A technique allowing for a perturbative treatment of nonlocal corrections to the single-site dy- 
namical mean-field theory (DMFT) in finite dimensions is developed. It is based on the observation 
that in the case of strong electron correlation the one-electron Green's function is strongly spatially 
damped so that its intersite matrix elements may be considered as small perturbations. Because 
the non-local corrections are at least quadratic in these matrix elements, DMFT in such cases may 
be a very accurate approximation in dimensions d = 1-3. This observation provides a rigorous 
justification for the application of DMFT to physical systems. Furthermore, the technique allows 
for a systematic evaluation of the nonlocal corrections. This is illustrated with two explicit exam- 
ples. First we calculate the magnetic short range order parameter for nearest neighbor spins in the 
insulating phase of the half filled Hubbard model on the square lattice which exhibits an excellent 
agreement with the results of a recent cluster approach. As a second example we study the lowest 
order correction to the DMFT self-energy and its influence on the local density of states. 

PACS numbers: 71.10.Fd, 75.10.-b 
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I. INTR OD UCTION 



The description of systems of strongly correlated electrons provided by the dynamical 
mean field theory (DMFT) has significantly improved our understanding of this class of 
materials. ^"^ The major limitation of DMFT lies in its single-site nature,^'^ which makes 
it unable to account for spatial correlations in finite-dimensional systems. To remedy this 
deficiency, cluster generalizations of DMFT have been devised (sec discussions on this sub- 
ject in Refs. 5 and 6). However, in order to recover the thermodynamic limit, clusters of 
sufficiently large size should be used.^ This strongly enhances the numerical effort needed 
to solve the corresponding DMFT equations, because the number of quantum variables to 
be simulated grows proportionally to the number of atoms in the cluster. To alleviate this 
difficulty, less costly but less accurate numerical approaches have been proposed. ^"^ In our 
opinion, this is not a viable solution, because the reduction of errors brought about by the 
larger number of atoms in the cluster can be overcompensated by the cruder treatment of 
the quantum effects. 

The aim of the present paper is to propose a perturbative approach accounting for nonlo- 
cal corrections to DMFT, with the latter being considered as the zeroth order approximation. 
The approach derives from the general method proposed in Ref. 10. This so-called gamma- 
expansion method (GEM) proved to be successful, inter alia, in the treatment of the local 
atomic correlations or short range order in the electronic theory of alloys. Here we 
apply it to the Hubbard model on a square lattice and in the insulating phase, which is well 
described by the alloy analogy (or Hubbard HI approximation) and therefore should be 
amenable to the same technique. This will be illustrated below with an explicit calculation 
of the magnetic short range order (SRO) for this model. The non-local magnetic correlations 
are of paramount importance for the theories of high temperature superconductivity based 
on the spin fluctuation mechanism, and their enhancement should therefore be of interest 
in these theories. In particular, both the Neel and the superconducting transition tempera- 
tures should monotonically depend on the strength of the short range magnetic correlations. 
In this respect the observation made in Ref. 9 that different cluster theories give short range 
correlations differing roughly by a factor of two provides an additional reason for developing 
non-cluster based treatments of non-local correlations. 

The proposed expansion is based in part on the widely accepted opinion that the single- 
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site DMFT which is exact in infinite dimensions provides us with a picture of the strong 
correlation phenomena which is quahtatively vahd also in physical dimensions d =1-3. In 
other words, it is the single-site dynamics which contain the essence of the problem while 
intersite correlations play lesser role and so may be considered as correction terms. 

An approximate non-local theory in finite dimensions, in which the single-site DMFT 
is used as the source of information on the local correlations has recently been proposed 
by Kusunose^^ and Toschi et al.^^. In essence, the method consists in first finding local 
irreducible vertex functions from the full DMFT correlators and then using these in the 
Bethe-Salpeter^^ or parquet^^ equations coupled with Dyson-like equations for the electron 
self-energy with nonlocal electron propagators in order to find the nonlocal self-energy. From 
our point of view the main deficiency of such approaches is the lack of consistency: The 
electron propagators in a skeleton expansion are the same in side-pieces and inside the 
crosspieces of a ladder diagram. Taking into account that for large values of the coupfing 
constant there is no reasons for the smallness of the terms in the ladder sums, in a strict 
mathematical sense the series are divergent. So the neglect of non-local contributions in 
some parts of the diagrams not only is unfounded but may also lead to serious numerical 
errors. In Ref. 15 a reference to the spirit of the DMFT in infinite dimensions was invoked 
where the irreducible vertices inside the ladder series can be chosen to be local without loss 
of accuracy.^'^^ This, however, is completely due to specifics of the d — oo case, while in 
finite dimensions the leading nonlocal corrections are are proportional to 0{l/d^^^y^ and 
thus are quite large in physical dimensions d =1-3. In our approach, however, the smallness 
of the non-local contributions can be explicitely established a priori. 

II. THE MODEL 

We are going to study the Hubbard model defined by the Hamiltonian 



where the number operators Uis for electrons with the spin projection s =t) i S'^'S expressed 
through the creation and annihilation operators cj^ and cLc,,; the intersite 

matrix elements tij connect nearest neighbor sites and coincide with the hopping integral t 
and the on-site matrix elements ta are set equal to the chemical potential /i. 




(1) 
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In the functional integral formalism (see, e.g., Ref. 18) the generating functional of the 
Green's functions is [note the absence of the caret on the Grassmann variables below in 
contrast to the operators in Eq. (1)]: 



Z[A, A+] = y" r>V exp I - j'^ 



dr 



+ H{r) + Y,At{Tmr) + Y.^t{r)A{r) 



(2) 



[cf. Eqs. (29) and (30) of Ref. 2]. Here the Grassmann variables A, A'^ are the spinor source 
fields which are conjugate to the physical spinor fields corresponding to the electrons 



and 



[ciT(^)>Cii(r)], 



(3a) 



(3b) 



where r is the thermodynamic imaginary "time" , and (3 = 1 jhsT., where T is the tempera- 
ture and kB the Boltzmann constant. Here and below we shall omit the constants originating 
from normalizations of functional integrals, because in the present paper we are interested 
only in correlation functions. The latter are obtained as functional derivatives of the gen- 
erating functional (2) at zero value of the functional arguments, divided by Z[0, 0], so the 
normalization is irrelevant. 

If the system is in its normal (i. e., non-superconducting) state and, besides, is paramag- 
netic, (the case we consider in the present paper), the exact one-electron Green's function 
Gijsij — t') as well as the self-energy T.ijs(T — r') arc independent of the spin direction 
* =T) i) so we omit this subscript for brevity. G and E are connected via the usual relation 
which we symbolically write as 

G^{dr + i-j:)-\ (4) 

The standard diagrammatic analysis (sec, c. g., Ref. 18) allows one to represent the 
generating functional of the Green's functions (2) as: 



Z[A,A+] = exp{A+GA)R[GA, A+G]. 



(5) 



To simplify notation, here and below the summations over the site and spin indices as well as 
integrations over the imaginary "time" are implicitly assumed in all products of quantities 



which depend on those mentioned. G in Eq. (5) is the exact Green's function defined in Eq. 
(4). The exponent on the right hand side of Eq. (5) corresponds to the free propagation of 
electrons while the functional R describes their mutual scattering and, diagrammatically, is 
represented by the matrix elements of the ^'-matrix with GA "tails" attached. 

III. THE METHOD 

As is well known (see review article 2 and references therein), in the infinite dimensional 
case the DMFT equations are exact. In finite dimensions DMFT may be viewed as a 
single-site approximation similar to the correlated effective field theory of ferroelectrics^^ 
or the coherent potential approximation in the theory of disordered alloys. ^° The latter 
approximations were found to be quite successful, so the non-local corrections arc expected to 
be small. In Ref. 10 a method of explicit calculation of these corrections was proposed, based 
on the observation that in the self-consistent perturbation theory the Feynman diagrams are 
expressed through the exact correlation function (or propagator) G. Thus, if the propagator 
is such that its on-site (or site-diagonal) matrix element is much larger than its off-diagonal 
elements, then, firstly, the single-site DMFT-type approximation should be a good one and, 
secondly, the dominant correction to this approximation can be calculated explicitly by 
taking into account those off- diagonal matrix elements which are the next largest after the 
diagonal ones. 

It was found that an adequate formalism for the realization of the above intuitive idea 
is provided by the functional-differential quantum-theoretical formalism. ^° In this approach 
the expression for the generating functional of the S'-matrix R in Eq. (5) can be cast into 
the form^^ 



which would coincide with the Hori representation^^ in the case of E = 0. In the above 
formula we have used the freedom in the separation of the total action [the term in the 
exponent in Eq. (2)] into the "free" (bilinear) part and the interaction part by chosing the 
(unknown) exact inverse Green's function G^^ as the free part. Then according to Eq. (4) 
this should be compensated by adding the term -0+2-0 to the interaction part Hu. The 
validity of the above representation can be easily shown both formally^^ and through a 
diagrammatic analysis. 
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The equation for the unknown exact one-electron Green's function can be obtained from 
its definition via the double functional differentiation of the generating functional (2): 

G = Z[Q,Q]-^5^hiZ[A,A+]/5A5A+\A,A+=o- 
Substituting here the expression (5) we obtain the following equation: 

5^hiR[^l^,'il^+]/5'il^i{T)5'il^nT') =0. (7) 



A. The gamma expansion method (GEM) 
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The last equation is exact and completely general (i. e., independent of any approximation 
scheme). To solve it in the case of strong correlation (large U) we first note that in Ref. 23 
it was shown that at finite temperature the one electron Green's function is exponentially 
damped at large distances. This confirms the original suggestion of Ref. 10 about the spatial 
asymptotic behavior of the correlation functions as 

G,,-ocexp(-|i-j|/0=C'(7l^"''l), (8) 

where 7 = exp(— 1/^). Assuming that the correlation length ^ is small, the Feynman 
diagrams containing large distance Green's functions can be neglected and the solution can 
be approximated by the local contribution plus a few terms accounting for correlations with 
neighboring sites. ^°Prom Eq. (8) it is easy to understand why our approach has DMFT as 
the zeroth order approximation. DMFT is exact in infinite dimensions because in dimension 
d the largest off-diagonal matrix elements of the Green function are bounded by a constant 
of 0(1/0?^/^) and so vanish at d = 00.^'^ The condition (8) with small 7 presumes that the 
off-diagonal elements nearly vanish. Thus, DMFT should be a good approximation in such 
cases in any dimension. 

Because the spatial dependence of the self-energy can be expressed through its dependence 
on the Green's function, we may assume that both G and E, considered as matrices in the 
lattice site indices i.j, arc diagonally dominated, so their separation into diagonal and non- 
diagonal parts is the separation into the leading contribution and correction terms: 

G^gl + G, (9) 

where / — Sij and g = Gu. Similarly, 

E = (77 + E. (10) 
6 



It is easy to check that if in the expressions given in the Appendix A one neglects the tilded 
(i. e., intersite) matrix elements of the Green's functions and of the self-energy operator, 
then Eq. (A4) substituted into Eq. (7) yields the DMFT equation.^ 

IV. MAGNETIC SHORT-RANGE ORDER 

To illustrate the above general approach we consider the simple problem of the correlation 
between electron magnetic moments at nearest neighbor (NN) sites on the square lattice at 
half-filling in the insulating phase. 

Following Ref. 9, we measure the magnetic moment in Bohr magnetons so, for example, 
the local moment at site i is = rii^ — rii^ = c\f^Ci, where is the Pauli matrix. In the 
paramagnetic phase, instead of the {m^rrij) correlator used in Ref. 9 we may consider the 
equivalent correlator 

(mr(r)m+(0)) = 2K(T)m^^(0)) (11) 

which is easier to calculate because it contains only the susceptibility function with opposite 
spins xn while the former correlator requires also the equal spin susceptibility. This is easily 
seen from the expression 

m- = {mty = 2c|,q, (12) 

which follows from the usual definition of the Pauh matrices = T^zkir'^. 

We define the magnetic short range order (SRO) parameters {mirrij) (we omit the super- 
scripts because they can vary in different definitions) as the zero-frequency Fourier compo- 
nent of the correlations functions (11) with i j which is equivalent to their average over 
the interval < r < /3. Because the statistical average of the product of the operators (12) 
contains at least two Green's functions (see Fig. 1), the lowest order contribution to the 
SRO parameter is 0(7^) [see Eq. (8)]. 

A formal expression for {mirrij) can be obtained from Eq. (5) by taking the appropriate 
fourth order functional derivative with respect to the source variables A and . The 
derivative of R[GA, A'^G] should further be expanded to second order in 7. The diagrams 
containing the contributions up to 0(7^) are shown in Fig. 1. The first diagram comes 
from the first factor in Eq. (5) corresponding to the disconnected contributions. The second 
diagram comes from the zeroth order (DMFT) term of the expansion of R and the last 
diagram corresponds to the 0(7^) term in the expansion of R (the two four-point single-site 
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FIG. 1: Diagrams contributing into the magnetic short range order to the lowest 0(7^) order of 
the small parameter 7. i and j are two nearest neighbor sites. The arrowed lines correspond to the 
fully renormalized one-electron Green's function G and the shadowed rectangles correspond to the 
fully renormalized amputated single-site vertex x^i (see Appendix A for the general definition) . 

vertices connected by two intersite Green's function lines). Here it is important to note 
that the diagrams of Fig. 1 should not be confused with the conventional decomposition 
of two-particle Green's functions into ladder-type series with the shaded rectangles being 
the irreducible four-point vertices as is the case with the diagrams shown in Fig. 8(b) of 
Ref. 2. Despite of their formal similarity, the meaning of the diagrams in the two cases is 
very different. While in the above-mentioned Fig. 8(b) the shaded rectangles in the case of 
physical dimensionalities would correspond to irreducible vertex functions with their spatial 
dependence fully taken into account, the shaded elements in our Fig. 1 are fuUi. e., in general 
reducible two-particle Green's functions but with their spatial dependence being only local, 
i. e., restricted to one site {k or /) even in physically relevant finite-dimensional case. The 
ladder-like structure of the diagrams shown in Fig. 1 is purely accidental and will not persist 
in higher orders. 

A. The Hubbard III approximation^^ 

To explicitly calculate the above diagrams we need expressions for their elements: the 
Green's functions and the four-point vertices. As can be seen from Fig. 2, the lowest order 
non-local correction to the DMFT sclf-cncrgy a is of order 0(7'^), so in calculations up to 
0(7^) the exact self-energy S may be approximated by a. 

For the Hubbard model at half-filling in its insulating phase the Fermi liquid properties 
are not important and the simple Hubbard HI approximation can be used.^^ Its detailed 
discussion, the connection with disordered alloys and the coherent potential approximation 
(CPA) as well as further references may be found in Refs. 20 and 2. Here we only mention 
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FIG. 2: The leading nonlocal correction to the DMFT self-energy. The meaning of the elements 
of the diagram is the same as in the previous figure. 



the observation made in the latter reference (see sect. VII. C. 2) that the Hubbard III ap- 
proximation can be obtained as a simple generalization of the atomic case as follows. Let 
us for generality consider the exact atomic Green function at an arbitrary occupancy n^^ 

I — n/2 n/2 

Qatiiuj) = ■ h -. ■ = Dat[iuJ + /X - (Tat{iuj)]. (13) 

lu) + II lu) + II — U 

If in this equation we replace the inverse unperturbed atomic Green's function by its DMFT 
counterpart as 

f Oat (^^) = 9at (^^) + ^at {i^) 

where, as usual, ^ 

Qo^iiu) ^ g~^{iu;) + a{iu;), (14) 

the Hubbard HI approximation 

/• N l-n/2 n/2 , , 

9sM = ^T-^ + \ (15) 

is obtained. 

Because the Mott transition takes place only at half- filling corresponding to n = 1, the 
insulating Hubbard HI solution is vahd only in this case. Furthermore, because at half filhng 
the system is explicitly particle-hole symmetric, it is convenient to switch to an explicitly 
particle-hole symmetric formalism. This is achieved through the replacement (see Ref. 2 Ch. 
VII B) 

^o"' ^ ^o"' + U/2. 

In this case Eq. (15) takes a symmetric form which we will use in the rest of the paper 

1/2 1/2 



9s{ii^) 



g^\iuj) + U/2 g^\iuj)-U/2 
= {l/2)g,,,{iuj) + {l/2)g,,_,{iuj). (16) 



Here, in the second line, we have introduced the partial single-site Green's functions 
corresponding to the propagation of electrons with their spins parallel or opposite to the 
dominant spin on the site under consideration (see Ref. 20). These Green's functions will 

be necessary in the next subsection. 
Eq. (14) changes to 

g^\iuj) = g-\iuj) - U/2 + a{iuj). (17) 

Eqs. (16) and (17) constitute a set to find g{iuj) and a{iuj), provided the bare single site 
Green's function is known. ^ The latter is given by the site-diagonal matrix elements of 
the lattice Green's function and is usually calculated as a sum over quasimomenta, which 
coincides with the Watson integral in the continuum limit^°: 

D{iuj) = ^^y-. ^— , (18) 

^ ' N ^ lu- e(k) ^ ' 

where e(k) is the lattice Fourier transform of the matrix —tij with = 0. With the use of 
this function the DMFT Green's function is obtained as g{iuj) — D[iuj + /i — cr(ia;)].^ 

Similarly, the matrix element Gi of the one-electron Green's function connecting NN sites 
can be calculated as Gi — Di[iu! + /i — a{iu!)], where, for isotropic matrices tij connecting 
only NN sites: 

Di{iu) = [1 - iujD{iuj)]/qt, (19) 

and q is the coordination number {q — 4: for the square lattice). 

We found it convenient to chose among various forms of the Hubbard HI (or CPA) 
equations^° the following one: 

a{iu) = Un/2 + {U/2fgo{iuj). (20) 

At weak couphng, when the first (Hartree) term in Eq. (20) dominates, the DMFT ap- 
proaches the perturbative solution. 

The solution for n = 1 was obtained by iterating Eq. (20) together with Eq. (17), and 
is shown in Fig. 3 for a particular set of parameters. As we see, the imaginary part of 
iuJn — o'{iUn) is quitc large, so the damping with distance of Gij in Eq. (8) is strong^^, hence 
7 is really a small quantity. For example, with the parameters given in Fig. 3 

max\Gi{iujn) \ ~ 2.5 • 10"V\ (21) 
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FIG. 3: The DMFT self-energy calculated in the Hubbard III approximation (present work) com- 
pared to the quantum Monte Carlo (QMC) simulations of Ref. 9 [un = 7rkBT{2n + 1) is the 
Matsubara frequency]. 

the correction terms being at least of the second order in Gi. This means that in the case 
under consideration DMFT is a very good approximation. 

In the next section we will see that in the case of strong correlation near half fiUing there 
exists a mechanism of enhancement of the loop diagrams at low temperatures. At such 
temperatures the 0(7^) non-local corrections corresponding to the diagram in Fig. 2 need 
to be taken into account. We will return to this question after explaining the mechanism of 
the enhancement. 

B. The Curie behavior of the magnetic SRO at low temperatures 

From Eq. (21) it follows that deviations from DMFT are small unless there exists some 
mechanism of enhancement of the correction terms. In this subsection we are going to show 
that in the half filled strongly correlated Hubbard model such a mechanism does exist and 
is connected with a peculiar time dependence of the four-point vertex functions entering the 
diagrams in Fig. 1. The term exhibiting this dependence is the first one in the expression 
for the atomic vertex function Eq. (Bl) and, at low temperature 

~ ^Su;i^Ju^^,j^Ag{iuj3)Ag{iuj4), (22) 
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where the superscript 'c' stands for 'connected', and 



Ag{iuj) = gs,s{iuj) - gs-s{ioj) = -Ugo{iu;)g{iu). (23) 

All other contributions to do not lead to an enhancement and can therefore safely be 
neglected. By analogy with Eqs. (13) — (16) above we may generalize the expression for the 
vertex function from the atomic to the band case as in the derivation of the Hubbard III 
approximation in the previous subsection, i. e., by replacing the atomic versions of gs,±s by 
their band counterparts. This can be justified more rigorously by a calculation of the vertex 
function in the Hubbard 111 approximation.^^ 

As was pointed out in Ref. 26, the singular character of the four-point vertex function 
is most transparent when written in the form of its contribution to the effective interaction 
functional (A7): 

where, according to Eqs. (22), (23), and (A8) 

(f){iu;) = Ugo{iuj)/g{iuj) (25) 

[g^^ appears due to the four amputations in Eq. (A8)]. It is easy to see that the above 
vertex function is a product of two independent bilinear terms with the symmetry of the 
S~ and 5"+ spin operators. When transformed to the (imaginary) time the interaction (24) 
takes the form 

^ I dTS-{T) J dr'S^ir'), (26) 

where the spin variables are the Fourier transforms of the two factors in Eq. (24) . Obviously, 
the above interaction can be exactly decoupled with the use of the Hubbard- Stratonovich 
transformation with a time- independent ordinary (i. e., non-functional) variable. This obser- 
vation provides a rigorous justification of the theories based on the static approximation. ^^"^^ 
We stress that, while in the theories of Rcfs. 27-29 the static approximation was introduced 
through an ad hoc modification of the initial Hamiltonian, in the present approach it was 
obtained as an approximate (Hubbard HI) solution to the Hubbard model, based on the 
observation made in Ref. 2 that, in the case of strong coupling and half filling, it accurately 
reproduces the results obtained in a fuU-fiedged DMFT treatment. In particular, our "by 
analogy" derivation of the vertex function can be confirmed by the observation that the 
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on-site spin-spin correlation function "develops long-term memory signalling the formation 
of a local moment" made on the basis of numerical simulations [see Ch. VII. G.l of Ref. 2 
and Eq. (B4) in Appendix B of the present paper] . 

Thus, because the static field has a zero-frequency Fourier transform we see that at 
strong coupling the half-filled Hubbard model develops a zero-energy scale. Below we will 
show that because of this the susceptibility may develop a Curie-like behavior on an energy 
scale much smaller than the electron bandwidth. The temperature dependence ~ 1/T will 
enhance the contribution of the corresponding diagram at low temperature, — thus providing 
a mechanism for the enhancement of the small 0(7^) contributions mentioned above. 

Formally this can be shown as follows. First, to the order 0(7^) we may neglect the 
difference between V and V in Eq. (A2) because according to Eq. (A3) their difference 
depends on the 0(7^) quantity E. Thus, the 0(7^) contribution to R in the spin fluctuation 
channel is obtained from Eq. (A2) as the 0{&) term in the expansion of the first exponential. 
The contribution (24) to V will then lead to the the following form for the dominant term 
in R 

SR = -a{T) c^^{(jj)(p{i(jj)cji{uj) ^ c+ (a;')(/)(iw')cfct(a;'), (27) 

where j ^-nd k denote nearest neighbor sites and 

a{T) = 5][Gi(iu;)^o(ic^)/^(ic^)]'. (28) 

Here the minus sign appears because of the closed Fermion line made by the two nearest 
neighbor Green's functions G^. The remaining factors come from the product of two effective 
interactions 4>{iu!). We note that the summation over the Matsubara frequencies in Eq. (28) 
is not supplemented by a f3~^ factor (this will be explained in more detail below). Thus, 
because the icu-dependent summand consists of smooth functions and the number of the 
Matsubara frequencies a;„ scales as the inverse temperature, as the temperature is lowered 
the sum will diverge as 1/T. The importance of this observation lies in the fact that this 
term leads to the magnetic susceptibility satisfying the Curie law which we are going to 
show with explicit calculations below. 

Among the diagrams of Fig. 1 it is the third one which exhibits the Curie behavior. This 
can be seen as follows. Because our magnetization operators Eq. (12) are dimensionless, 
their correlator (11) or (30) also has no dimension. The elements of the diagrams in Fig. 1, 



13 



on the other hand, when Fourier transformed with respect to the imaginary time variables 
have the following dimensionalities in units of energy: 

= and = [x?J[G]-^ = 2 (29) 

[see Eq. (Bl)]. Thus, in the first diagram in Fig. 1 these elements contribute the factor 
GG which has dimension -2. Because the contribution should be dimensionless, this means 
that the diagram additionally contains the compensating factor {ksTY with the dimension 
2. According to the diagrammatic rules corresponding to the vertex function (24), in the 
particle-hole channel the number of summation over the internal variables of a diagram is 
equal to the number of loops in this diagram. The summation over the Matsubara frequencies 
scales as ~ 1/T for each loop as explained above. Thus, in the case of the first diagram 
we have one loop summation of 0(1/T) and the factor O(T^) so the result is bounded as 
T ^ 0. The second diagram has the combination GGx'^GG of dimension -2, hence the same 
0{T'^) additional factor but this time two loop summations bringing the 0(1/T^) behavior 
making the diagram to be of order unity. Finally, a similar analysis in the case of the last 
diagram in Fig. 1 shows that it behaves as 0(1/T) because of its three summations over the 
loop variables. It is this latter diagram which are going to compute explicitly below while 
neglecting the other two because of their smallness. 
The contribution corresponding to this diagram is 

{m^ rrij) = (5 ^ j dT{m^ (r)mj^(O)) 
Jo 

— mimja{T), (30) 
where the on-site magnetizations are 

^iij) = -P'' = up-' 9{^^)9o{^to). (31) 

To understand the origin of the contributions to Eq. (30) we first note that the two lines 
connecting sites k and / on the diagram come from the intersite correction to the four-point 
vertex in Eq. (A2) where to the order 0(7^) we may restrict G to the NN sites. This 
explains the appearance of G^ in Eq. (28). Thus, the smallness of the diagram in Fig. 2 
is already exhausted by this contribution, so the leftmost and the rightmost propagators 
may be restricted to on-site Green's functions g. In other words, to the order 0(7^) we 
have i — k and j — I. The factors g'^ at sites i and j multiply the two external legs of the 
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vertex (22) which, combined with Eq. (25), leads to expressions of the type of Eq. (23). The 
summation over frequencies then yields the two factors of the type of Eq. (31). The on-site 
magnetizations mi(j) (31) calculated with the parameters of Fig. 3 are practically saturated 
to their maximum value 1 (in Bohr magnetons). 

In the two internal lines the amputation factors are not compensated by the propagators 
because they are now intersite ones (Gi), so the two vertices (24) introduce the two fac- 
tors Ugo/g. The coefficient (1/2) x (1/2) coming from the product of two vertices (24) is 
compensated by the coefficient 2x2 coming from Eq. (12). 

The magnetic SRO calculated according to the above formula shows excellent agreement 
with the cluster calculations of Rcf. 9 (see Fig. 4). Because our technique works directly 
with the thermodynamic limit, this agreement confirms the conclusion of the above reference 
that their Active impurity (FI) method is a more reliable cluster approximation than the 
dynamical cluster approximation which yields a SRO parameter approximately two time 
larger than in Fig. 4 (see Fig. 14 in Ref. 9). 

Because of the 0{1/T) enhancement, our perturbative technique becomes less rehable 
at very low temperature, where a resummation technique should be used. For example, by 
extending the ladder sum to infinity with the intermediate Green's functions restricted to 
Gi as above, we would have obtained the usual mean field value for the Neel temperature 
Tn — qj. The value of calculated in this way would match that shown in Fig. 10 of Ref. 
9. Such an approach, however, is not in the spirit of the GEM. In fact, within the GEM the 
phase transition temperatures can usually be calculated with much better accuracy than in 
the mean field approximation. For example, according to Eq. (13) of Ref. 10 the value of 
gJ/Tjv in the zeroth order approximation of GEM is given by the Watson integral for the 
corresponding lattice. In a 2-dimensional system this integral would diverge (see, e. g., Ref. 
30), — thus giving the correct value = in this case. The possibility to develop such an 
approach for phase transitions in the Hubbard model is currently being investigated. 

A Gurie-like behavior is characteristic for localized spins but is difficult to understand 
in the case of itinerant spins of band electrons. The latter situation may be described by 
the Hubbard model in the weak coupling limit U ^ 0. Because in this case conventional 
perturbation theory applies, all calculations can be done in the finite bandwidth case. But 
it is more instructive to see how the perturbative limit is attained in the formalism of the 
GEM. To this end we first note that the approximation (22) is not appropriate in the weak 
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FIG. 4: Magnetic short range order parameter for the NN sites i and j calculated to the lowest 
nontrivial order of the gamma expansion method (GEM) in comparison with the cluster technique 
of Ref. 6 [called by the authors the fictive impurity (FI) method] and with the leading term of the 
high temperature expansion of the superexchange interaction. ^'^^ The square dots correspond to 
data calculated by two different methods in Ref. 9. 

coupling case because it is of second order in U . Going back to the defining equation (Bl), 
we find that the first order term has the small- C/ limit 



where we presented the amputated version of the connected correlator because in this case 
the atomic correlator is the same as in the band case. As we see, now the vertex functions 
in the diagrams of Fig. 1 contain an additional factor T, so that the mechanism of the 
low-temperature enhancement leading to the Curie law is not operative any more. Because 
of this and because of the smallness of U it is the first diagram in Fig. 1 which is the main 
contribution with the second diagram giving the first order correction. Because the first 
order Hartree term in the Green functions is exact and the vertex in Eq. (32) is also exact, 
the spin susceptibility calculated within the GEM will be accurate to first order in U at 
weak coupling. 

C. Non-local corrections to the self-energy 

The mechanism at the origin of the enhancement in the loop summations in the half- 
filled Hubbard model will also be operative in the diagram for the lowest-order non-local 




(32) 
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correction to the self-energy shown in Fig. 2. Indeed, it is easy to see that (i) the contribution 
to E coming from the term (27) can be obtained by simply contracting Cis and at nearest 
neighbor sites and (ii) by noting that the contribution thus obtained should be multiplied 

by 3/2 because the term (27) coming from the product (26) constitutes only 2/3 of the full 
scalar product S ■ S = S^S^ + S^'S^ constituting the full rotationally invariant singular 
interaction^^ leading to the first term in Eq. (Bl). Thus 

Ei(z^) = ^[UgoM/gM]MT)Gi{tuj). (33) 

In connection with this expression some comments are in order. Prom Eq. (28) it follows 
that the nonlocal self-energy in Eq. (33) is of fourth order in U . Such a term was shown to 
be the lowest order (in U) non-local contribution to the self-energy in the Falicov-Kimball 
(F-K) model in Ref. 33. Our approximation (22) essentially describes this case because it is 
based on the local (atomic) approximation to the band electrons — a feature shared by the 
/-electrons in the F-K model.^^ In the Hubbard model, the last, linear in U term in Eq. 
(Bl), dominates in the weak coupling limit. Its insertion into the diagram of Fig. 2 gives 
the conventional second order contribution into S. 

D. Herglotz properties of spectral functions 

Our next comment concerns the fundamental requirement that the imaginary part of the 
self-energy be negative semi-definite in order to guarantee the positivity of the density of 
states 

p{E,k) = --QJ:{E,k)\G{E,k)\' (34) 

TT 

of the quasiparticles with quasimomentum k and energy E. The quantities in the right 
hand side of this equation should be obtained by analytic continuation from the discrete 
Matsubara frequencies in the upper part of the complex energy plane z to the real axis E as 

z = E + ie\^^^^. 

This continuation is in general a non-trivial procedure because the quantities to be contin- 
ued are usually known only in the form of the data of quantum Monte Carlo simulations. 
Continuation of such data to the real axis is an ill-posed problem the discussion of which as 
well as the pertinent bibliography may be found in Ref. 2. We only would like to add that 
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the continuation of a perturbation expansion (like GEM) poses additional problems which 
are explained in Appendix C. 

Because of the simplicity of the Hubbard III approximation, however, it is possible to 
perform the above analytic continuation explicitly. This will be done below to illustrate 
the Herglotz properties of the self-energy and the total density of states with the correction 
term given by Eq. (33). To this end we first derive the expression for the corrected on-site 
Green's function. For a hopping matrix connecting only nearest-neighbor sites, we get from 
Eqs. (4) and (18)] 

g{iuj) = Gii{iuj) = D . (35) 

\ — l^i{iuo)/t \ l — l^i[iuj)/t J 

From this expression it is easy to see that irrespective of the precise form of the correction 
term, the total density of states (DOS) 



p[E) = --%(^ = ^ + ie) 

TT 



(36) 

e->-0 

is conserved because the asymptotics of g{z) at large z are the same as those oi D{z): ~ 1/ z 
(sec Appendix C for further details). However, from the point of view of GEM, Eq. (35) is 
valid only to lowest order in Si. So it is more prudent (and also more convenient for further 
analysis) to separate different contributions as 

g{iuj) ~ D[iu — aiiuj) + /i] 

+ [Mi^)/t\ [D{z) + zD'{z)l^,^_^^,^^^^ , (37) 

where the second line represents the lowest order correction to DMFT (the first line). 

Eqs. (20), (33), (35), and (37) were continued to the real axis and solved with the model 
parameters used in the calculations above. The results did not show any violation of the 
Herglotz properties. Therefore, to study the anomahes discussed in Appendix G, a smaller 
value olU = 3.2t was chosen, which corresponds to 80% of the critical value Um ~ 4:t for 
the Mott transition in the Hubbard 111 approximation, and leads to a more than 30 times 
larger NN matrix element (21) . 

Before proceeding with explicit calculations we would like to remind that our approach 
is based on the asymptotic behavior of the imaginary-time Green's functions at finite 
temperature. The spatial attenuation (8) is defined by the Green's function behavior at 
the discrete imaginary energies (the Matsubara frequencies), being more efficient at larger 
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temperatures.^^ Therefore, there is no guarantee that the correction terms will remain small 
when continued to the real energy axis. Explicit calculations below, however, show that the 
GEM corrections remain generally small even on the real axis. 

The GEM-corrected self-energy in energy and wave vector space can be written as [cf. 
Eq. (10)] 

E(^,k) a{z) + 2(cos/i;a; + cos ky)a'E[{z), (38) 

where the prime means that the self-energy (33) has been divided by the SRO parameter 
a{T), which is the natural expansion parameter in the present case. Since a{T) is negative, 
Eq. (C3) takes the form 

2( cos kx + cos ky)J2'JE)'\~^ 

max ^ , , ^ . 39 

aE'i(E,k)>o a{E) J 

An explicit calculation shows that for the case under consideration the imaginary part of 

E(2;, k) is non-negative only when a > —0.45. We do not consider this restriction as a severe 

limitation of the GEM because, in such a situation, the spins in the system are spatially 

so strongly correlated, that the behavior of the electrons changes qualitatively and more 
sophisticated techniques should be developed for calculating the Green's functions (see, e. 
g., Ref. 36). 

The results of our calculations are shown in Fig. 5. The Hubbard sub-bands a,t U — 
±OAUm are not yet separated. Prom the shape of the derivative of the DOS with respect 
to the SRO parameter (lower panel), we expect a transfer of spectral weight from the sides 
of the subbands to their centers as the temperature is lowered and a{T) becomes more 
negative. The physics of this transfer is the same as in the Falicov-Kimball (F-K) model of 
Ref. 35. When the local order grows, the atoms became increasingly surrounded by sites 
with essentially different on-site potentials. In the F-K case these are the atoms filled by the 
/-electrons while in the Hubbard case the surrounding atoms host the electrons of opposite 
spins. This hampers intersite atomic hopping and causes the subband narrowing. In the 
F-K model case this can be seen from Fig. 7 of Ref. 35. The qualitative difference with 
our Fig. 5 is only due to (i) the narrow features at the centers of the subbands because 
of the atomic-like spectrum of the /-electrons which are absent in the Hubbard model and 
(ii) stronger temperature dependence of DOS near the middle of the band because in the 
2D F-K model a charge transfer gap opens in the the spectrum at finite temperature. This 
phase transition at finite temperature is possible because the order parameter of the charge 
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FIG. 5: Upper panel: Density of states (DOS) of the half-filled Hubbard model on the square 
lattice with U = 3.2t calculated for two values of the magnetic short range order parameter a. 
The latter was calculated with the use of Eq. (28) for the two temperature shown in the figure. 
The DOS was calculated via the first order expansion Eq. (37). Lower panel derivative of the DOS 
with respect to the short range order parameter. All quantities are measured in energy or inverse 
energy units, with the the hopping parameter t set equal to unity. 

density wave is Ising-like, while in the 2D Hubbard model the ordering is forbidden due to 
the continuous (rotational) symmetry of the spin variables. 

Finally, it is worthwhile to note that the enhancement of the intersite correlations dis- 
cussed in this work for "normal" contractions also takes place for the anomalous means 
(cjCj) in which case it leads to a BCS-type superconducting gap equation with the antifer- 
romagnetic intersite correlations playing the role of the pairing interaction (see the review 
paper Ref. 14). This problem will be considered in more detail in a future publication. 

V. CONCLUSION 

As our elementary calculations show, in some cases the GEM allows for a reliable calcu- 
lation of non-local corrections to DMFT. In the case of the Hubbard model at half filling the 
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corrections are strongly enhanced by spin fluctuations. We expect that a similar (though 
presumably weaker) enhancement will also take place in the case of small doping. Fur- 
thermore, from the discussion above, it is seen that the enhancement occurs in the leading 
non-local correction to the four-point vertex (last diagram in Fig. 1) and to the self-energy 
(particle- hole loops in Fig. 2). In both cases the enhancement may be of considerable im- 
portance for the spin fluctuation mechanisms of superconductivity formulated either via the 
Bethe-Salpeter equation for the bound states of the pairs or the BCS equation for the self- 
energy. ^^'^^'^^ This would require more accurate calculation of the four-point vertices with 
proper account of the Fermi liquid properties. This problem is currently being investigated 
in the framework of the approaches of Refs. 2,39,40, and 41. 
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APPENDIX A: RENORMALIZED VERTICES 

The relationship between the conventional path integral representation of partition func- 
tions and the Hori functional-differential representation can be formally established by first 
using the Gaussian functional integration formula to decouple the second order functional 
derivative into a first order derivative multiplied by a Grassmann field (say, rj) in a way 
completely analogous to the Hubbard-Stratonovich transformation and then applying the 
formula for the functional shift (Ref. 18, Ch. 1), 



valid for any functional / (and similarly for 77+, ■0+). Next with the sequence of transforma- 
tions explained in Refs. 10 and 25 one obtains: 



exp{r|^/^'4))f[^l)] = /[^ + r]] 



(Al) 



i.[^,^1=exp(A6^ 



) 



exp(v^+rV- 



(A2) 
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(A3) 



+ rj+{g-' + a)rj-Hj[rj,rj+]}y 



where 



and the spinors 77 and 77''" have the structure of Eqs. (3a) and (3b) 



7]{t) -- 

and 

(r) 



a|(r) 



(A4) 



(A5) 



(A6a) 



(A6b) 



In Eqs. (A2) and (A3) the summation over the site indices and the integration over the 
thermodynamic "time" r is imphcitly assumed in all products of fields and matrices. In 
Eqs. (A4)-(A6b) r] is not supphed with the site index because we consider the normal, 
paramagnetic, and translationally invariant case in which all sites are equivalent. In this case 
the matrices g and a are diagonal and proportional to the unit matrix because g] = gi = g 
and similarly for the self-energy. The formalism, however, can be easily generalized to a 
general broken symmetry case which would essentially amount to introducing site-dependent 
field ?7 and the the Green's function and self-energy matrices of general form.^°'^^ 
Thus, according to Eq. (A4) the effective on-site interaction vertices are given by 




'^n — 1 '?~n 



even n > 2 

^i{'^i)'4}i{^2) ■ ■ ■ X^('ri, T2, . . . , T„_i, r„) . . . A{rn-Mi{rn), 



(AT) 



where the "amputated" vertex is defined by 



Xn(^i> '^2, ■ ■ ■ , T„_i, T„) = / / . . . / g \ri-T[)g \t2 -tQ... 



(A8) 



The summation in Eq. (A7) is restricted to the values of n > 2 because according to the 
zeroth order of GEM (or DMFT — see below) approximation Eq. (7) with i — j the term 
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n — 2 turns to zero, so V[ipi, contains only true many-electron interactions. The angular 
brackets in Eq. (A8) denote the statistical averaging in Eq. (A4) with the effective action^ 

Seff^V^{9-' + ^)v-Hi[v,V^] (A9) 

with the superscript 'c' meaning that only the connected diagrams should be kept. It should 
be noted that because ipi and 77 in the above equations are of the spinor type, the vertices 
Xn s-re tensor quantities possessing in general 2^" components. In the symmetric case under 
consideration many of the components are equal zero and many among them are mutually 
equal. For example, the n = 2 four-point function is fully defined by two scalar vertices: 
Xll and of which in our calculation of the spin correlation function only the first one is 
needed. 

Thus, the effective single site vertices in Eq. (A7) are the connected [because of the loga- 
rithm in Eq. (A4)] amputated [because of the factors in Eq. (A8)] single-site correlation 
functions calculated with the effective action given by Eq. (A9). In Eqs. (A4), (A5), (A7), 
and (A8) all fields and matrices are assumed to be single-site. To simplify notation the site 
dependence has been omitted. It can be easily restored by supplying all quantities in these 
formulas with the subscript i. 

There is an infinite number of spatially local interactions in the GEM but, in contrast 
to what happens in conventional perturbation theory, the smallness is provided by the off- 
diagonal (in the lattice site indices) elements of the propagator G and of the self-energy E 
[see Eqs. (A2), (A3), and (A4) above]. So the order of a correction will be classified by 
the number and length of the propagators in the corresponding diagram. For example, if 
Eq. (7) is restricted to a single site, i. e., if we assume that i = j, then Eqs. (A4), (A5) 
with all tilded quantities set to zero together with Eq. (7) will reproduce the DMFT theory 
of Ref. 2. Because i — j means the zeroth order of the GEM [see Eq. (8)] the above means 
that DMFT is equivalent to the zeroth order of the GEM. 

APPENDIX B: ATOMIC FOUR-POINT CORRELATION FUNCTION 

In this appendix we present the connected part of the atomic spin susceptibility which 
was calculated exactly in the Appendix of Ref. 42 where this quantity was denoted as Fo-.-o-- 
In the present paper we stick to more common notation of Ref. 2. We consider the limit 
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U » ksT and neglect all terms of the type exp(— [7/2/^^7"). Besides, the particle-hole 
symmetry is introduced by setting the one-electron energy e to —U/2. 

= {ai{iu;i)ai{iu;2)ai{iu;3)ai{iu;4)y 



VP 



e 



iuJlTl +iui2T2—i(jJ3T3 —lUJiTi 



{a^{Ti)ai{T2)a^{T2,)ai{Ti)Y 



fc=i V^o 

-U I 5a;iu;4(^a;2a'3 + 

+^kBTUv{{iu}k})S^^+^^^^^+^^, (Bl) 



2 



where 



and 



1 



[zur - {u/2y 



v{{iu}k}) = [ 2 JJ(ia;fe + U/2) - U{iui + iu2 + C/) 
\ it=i 

X [(icua + t//2)(ia;4 + U/2) + (icui - [//2)(icu2 - U/2)] \ \[ f{iujk). 



k=l 

The atomic approximation may look to be too crude to use in the finite band width case and 
in general this is true. Nevertheless, it captures in a qualitatively correct way the important 
particular cases which we used in our study. The most important to us is the first term in 
Eq. (Bl) which we used in our calculations in the main text. If factor out the kinematic 
5-function responsible for the total energy conservation as 



2 



-^SuJlU14^U!2UJ3f (i'^s) f i'i'^4) — SuJl+UJ2,UJ3+UJ4. 



2 



^S^2-uJ3,ofii'^3)f{^'^4:) 



(B2) 



one can see that the dynamics of this term is defined by the zero-energy delta-function 
singularity in the particle-hole channel corresponding to the symmetry of the S"^ spin op- 
erators. A zero-energy excitation in the spin channel is not surprising in the atomic case 
where it represents simply the Goldstonc boson due to the rotational symmetry breaking 
in the ground state in the half-filled case. Because of the rotational symmetry the electron 
conserves the direction of its spin. 

But in the case of itinerant electrons the magnetic moments at individual sites are not 
conserved. Still, the DMFT simulations of Ref. 2 showed that in sufficiently strongly corre- 
lated half-filled Hubbard model the 4-point correlation function does develop the 5-function 
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singularity in the particle-hole channel. This, — in particular, — is reflected in the (imaginary) 
time independence of the spin-spin correlation function [see Eq. (252) in Ref. 2]. This can 
be seen as follows. According to Eq. (Bl) [see also Eqs. (11), (12), and (23)] 

X^oti-r) ^ {m-{r)mtm ^ p-' ^ e^(--^-3)r 

U)I...U>4 

x25^^^^5^^^^/^g{iuj^)/^g{iuji) ^ 2, (B3) 

where in the last line we used the saturation of the magnetic moments at strong coupling (see 
the main text). Thus, the term under discussion describes static local magnetic moments 
because their correlation function is independent of r. Because of this, the Fourier transform 
has the only component which is different from zero 

Xr.+ (i^ ^^)= f drxiotir) ^ W (B4) 

which is approximately twice the result obtain in the QMC simulations of Ref. 2 [see their 
Fig. 44 and Eq. (254)] as it should be in view of Eq. (11). 

It should be stressed that the above results are valid only in the insulating phase of the 
model. The insulator phase appears for U larger than some critical value of C/ — Uc and 
only at half- filling. It is precisely in this case that we apply the approximation (22). 

For U < C/c in the band case the 5-function singularity smears out and transforms itself in 
one among many others 0(f/^) contributions into v in Eq. (Bl). In the atomic approximation 
it will persist at all values of coupling because the atom is always in the insulator phase. But 
at weak coupling the singularity will be dumped by the coefficient U"^ while the dominant 
linear in U term given by Eq. (32) is exact also in the finite band width case. Thus, our 
theory interpolates between two limits where it agrees with known reliable approaches. 



APPENDIX C: ON PERTURBATIVE EXPANSION OF SPECTRAL DENSITIES 

Let A{E) be a function corresponding to some spectral density, which means that the 
function must be (i) positive and (ii) be normalized to some constant value: 

/oo 
dEA{E) = C. (CI) 
oo 

The above properties are typical of various spectral functions of the Green's function theory, 
such as the densities of states (34) and (36) with C = 1 or the imaginary part of the electron 
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self-energy multiplied by — I/tt with C — U'^n{l — n) in the paramagnetic case (see, e. g., 
Appendix A of Ref. 43). 

Condition (ii) is a consequence of the asymptotic behavior of the type C / z as z ^ og 
of the corresponding complex function of which the spectral function is the imaginary part 
just above the real axis, multiphed by — l/vr. 

Because the higher order terms of a perturbative expansion usually contain the products 
of a larger number of Green's function than the lower order terms, Eq. (CI) is usually satisfied 
in every order of approximation, provided the order is sufficiently large. (In the case of the 
self-energy, for example, the order should be larger or equal to two.) In particular, the site- 
nondiagonal matrix elements of the Green function which wc use in the gamma expansion 
are always quantities of order 0(1/2;), so the correction terms of any order cannot modify 
the 0{1/ z) asymptotic behavior of a function under consideration. 

Now let us consider some approximation A^^^ > to the function A{E) satisfying Eq. 
(CI). It is easy to see that any correction term of order m, \^A^'^\E) (where A is the 
expansion parameter) should satisfy the restriction 

/oo 
dEA^"'\E)^Q. (C2) 
-00 

From here it follows that A^'^\E) acquires both positive and negative values. Let us for 
definiteness assume that A is positive (the case A < is treated similarly) . Now if 

i^(™)(E)r"' 



A"* > 

then, as is easily seen, the corrected 



max 



A(-)<o A(o)(E) 



A{E) ^ A'^^\E) + y^A^'^^E) = A'^^^E) 



a^'^He) 

1 _l_ ^ ' 

^ A(^){E) 



(C3) 



(C4) 



will acquire negative values. 

Thus, in perturbative calculations the non-negativity of spectral functions can be ex- 
pected only for sufficiently small values of the expansion parameter. For larger values special 
care must be taken (e. g., partial resummation of infinite sequences of the diagrams) for the 
calculated spectra to be physically acceptable. In general it should be expected that in a 
series expansion any exactly known property can be guaranteed to be accurate only up to 
the size of the next order correction provided the series converges. According to Eq. (C4), 
the condition (C3) corresponds to a situation where the correction term becomes larger than 
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the main term, which signals a possible divergence of the series. In that case a partial sum 
of the series cannot be considered as an approximate representation of the function and 
usually cannot reproduce its general properties, such as the positive definiteness. 
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